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1.  Introduction 

We  are  concerned  with  the  identification  and  characterization  of  solution  errors 
in  spherically  symmetric  shock  interaction  problems.  This  issue  applies  to  the  study 
of  supernova  and  the  design  of  inertial  confinement  fusion  (ICF)  capsules.  In  the 
first  case,  theory  and  simulations  contain  a  number  of  uncertainties,  and  comparison 
to  observations  is  thus  not  definitive.  A  systematic  effort  to  remove  some  of  the 
uncertainties  associated  with  simulation  will  thus  be  a  useful  contribution.  In  the 
second  case  of  ICF  design,  concern  over  solution  accuracy  has  led  to  mandates  of 
formal  efforts  to  assure  solution  accuracy. 

In  previous  papers,  we  have  developed  a  general  approach  to  uncertainty  and 
numerical  solution  error  [5,  6],  and  we  have  analyzed  shock  interactions  in  a  planar 
geometry  [2,  4].  Here  we  focus  specifically  on  complications  which  result  from 
spherical  geometry.  In  brief,  these  are: 

(1)  The  solution  waves  are  not  of  constant  strength  between  wave  interactions, 
but  evolve  approximately  according  to  a  power  law  as  a  function  of  the 
radius. 

(2)  The  solution  is  not  spatially  constant  between  waves. 

(3)  If  the  solution  is  not  required  to  be  spherically  symmetric,  the  problem  of 
identifying  wave  structures  as  curves  or  surfaces  in  2D  or  3D  is  introduced. 

These  are  classical  problems  as  far  as  the  solutions  are  concerned,  but  the 
application  of  these  ideas  to  the  analysis  of  errors  in  the  solution  appears  to  be 
new.  The  radially  dependent  strength  of  spherical  waves  is  discussed  in  [8].  The 
spatial  variation  of  spherical  waves  is  contained  in  the  Guderley  solution  [7], 

The  problem  of  analysis  of  errors  in  numerical  solutions  is  of  course  central  to 
numerical  analysis.  Much  of  this  effort  is  motivated  by  other  concerns,  and  appears 
not  to  be  directly  applicable  to  the  problems  we  address. 

An  early  focus  of  numerical  error  modeling  was  round  off  errors.  For  the  hyper¬ 
bolic  systems  we  study,  modern  64-bit  processors  with  double  precision  arithmetic 
appear,  els  a  practical  matter,  not  to  be  sensitive  to  this  class  of  errors,  while  they 
are  difficult  to  analyze  theoretically.  A  more  common  approach  to  error  analy¬ 
sis  in  numerical  analysis  is  the  study  of  the  asymptotic  behavior  of  errors  under 
mesh  refinement.  This  is  a  useful  approach,  and  one  we  refer  to  in  the  case  of 
well  resolved  simulations.  However,  we  want  an  analysis  which  is  also  applicable 
to  the  pre-asymptotic  case  of  under  resolved  simulations  as  these  are  so  typical  of 
practical  studies  of  realistic  complex  physical  systems.  Moreover,  the  coefficients 
which  multiply  powers  of  Ax  in  the  asymptotic  expressions  cannot  be  determined 
theoretically.  A  third  main  theme  in  the  analysis  of  errors  is  the  use  of  a  posteriori 
error  estimators.  These  are  upper  bounds  on  the  error  in  the  solution,  based  on 
criteria  derived  from  the  (approximate)  solution  and  the  exact  equation  only.  Such 
a  study  comes  close  to  the  problems  we  address,  but  differs  in  a  few  respects.  We 
seek  to  characterize  the  error,  not  just  to  bound  it.  Moreover,  a  posteriori  methods 
are  most  fully  developed  and  justified  theoretically  for  elliptic  problems,  and  have 
only  a  partial  or  preliminary  development  for  the  shock  interaction  problems  we 
consider  here.  A  fourth  approach  to  errors  is  to  regard  them  as  due  to  input  uncer¬ 
tainties.  In  this  point  of  view,  uncertainty  analysis  is  a  mapping  of  input  random 
variables  to  output  random  variables. 
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In  contrast  to  the  first  three  of  the  above  methods,  we  analyze  the  errors  sta¬ 
tistically.  The  statistical  models  are  simple,  and  in  this  sense,  we  analyze  only 
the  central  portion  of  the  error  statistics,  not  their  tails.  We  use  a  Gaussian  error 
model,  and  thus  we  identify  the  mean  error  and  the  covariance.  The  mean  error  is 
a  systematic  error,  and  it  can  in  principle  be  used  to  modify  or  “post  process”  and 
error-correct  the  approximate  solution.  The  covariance  is  a  measure  of  the  error 
variability.  One  may  question  the  idea  that  numerical  errors  can  be  modeled  statis¬ 
tically  or  that  the  errors  are  variable,  when  each  simulation  is  totally  deterministic. 
This  philosophical  question  has  an  easy  answer:  determinism  lies  in  the  eye  of  the 
beholder.  In  other  words,  the  modeling  of  a  natural  phenomena  (tossing  a  coin, 
for  example)  as  a  probabilistic  or  deterministic  event  depends  on  the  level  of  detail 
included  or  excluded  from  the  model.  Thus  we  argue  only  why  it  is  convenient  or 
nearly  essential  to  omit  from  the  modeling  of  error  information  needed  to  make  the 
error  analysis  deterministic.  The  variable  (as  opposed  to  the  systematic)  part  of 
the  error  depends  on  “accidental”  features  of  the  numerical  model,  such  as  the  sub 
grid  location  of  waves  and  wave  interaction  points  relative  to  the  mesh  cell  edges. 
Clearly  a  deterministic  model  requiring  such  data  would  be  too  cumbersome  for 
use  in  practice,  and  thus  a  probabilistic  model  is  preferable.  Indeed,  the  essentially 
probabilistic  aspect  of  round  off  error  is  well  recognized.  Given  that  the  complex¬ 
ity  of  a  statistical  model  is  needed,  we  found  that  simple  linear  error  models  were 
sufficient  for  our  analysis  [4].  The  simple  reason  for  this  pleasant  turn  of  events 
is  that  the  error  is  similar  to  a  perturbation,  and  normally  a  small  perturbation. 
Thus  the  error  of  a  strongly  nonlinear  problem  still  has  a  useful  linear  expression, 
at  least  as  far  as  our  analysis  of  the  error  has  progressed. 

In  contrast  to  the  fourth  approach  to  uncertainty  quantification  above,  we  allow 
for  errors  generated  within  the  solution  pROcesses.  Thus  we  subsume  and  expand 
on  this  point  of  view. 

As  a  technical  introduction  to  the  paper,  we  study  errors  of  a  spherical  shock 
interaction  problem  on  uniform  radial  grid  of  100  and  500  cells  (with  errors  deter¬ 
mined  by  reference  to  a  2000  cell  calculation  referred  to  as  the  fine  grid).  We  use 
MUSCL  [1]  as  the  numerical  method;  for  the  comparison  of  tracked  to  untracked 
solutions  of  the  problem,  see  [3],  The  equation  of  state  is  a  7-law  gas  with  7  =  5/3. 
The  ensemble  of  200  initial  conditions  is  defined  by  a  Latin  hypercube  variation 
shock  and  contact  strength  by  ±10%  about  a  base  case  defined  (as  in  [4])  by  a 
contact  located  at  1.5  units  from  the  origin;  an  inward  moving  shock  located  2.25 
units  from  the  origin,  with  all  constant  states  between  waves.  The  initial  base  case 
shock  strength  is  M  =  32.7  and  the  initial  base  case  Atwood  number  for  contact  is 
0.82. 


2.  The  Statistical  Numerical  Riemann  Problem 

We  study  statistical  numerical  Riemann  problems  (SNRP)  in  spherical  geome¬ 
try.  The  SNRP  introduces  errors  (modeled  as  random)  in  addition  to  propagating 
errors  or  uncertainty  from  input  to  output.  The  waves  in  the  SNRP  have  a  finite 
width  and  the  solution  algorithm  in  the  SNRP  has  only  finite  accuracy.  Because  of 
the  possible  finite  width  of  the  input  waves,  the  problem  and  its  solution  are  not 
strictly  scale  invariant.  Moreover  scale  invariance  is  lost  through  the  length  scale 
introduced  by  the  radius  at  the  time  of  interaction. 
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We  consider  a  generalization  of  the  spherically  symmetric  Riemann  problem.  A 
statistical  distribution  of  numerical  incoming  waves  and  starting  states  determines 
the  SNRP.  Its  solution  gives  the  output  waves,  each  of  which  generate  the  same 
type  of  data.  Thus  we  define  the  SNRP  as  a  statistical  (non-deterministic)  mapping 
from  a  statistical  input  wave  description  to  a  statistical  output  wave  description. 

The  statistics  of  the  SNRP  mapping  function  arise  from  grid  errors,  and  from 
the  random  placement  of  a  traveling  wave  relative  to  the  centers  of  the  finite  dif¬ 
ference  lattice.  Our  objective  in  this  section  is  to  build  up  a  library  of  statistical 
input-output  relations  that  will  include  all  Riemann  problems  to  be  encountered  in 
Sec.  3.  This  library  will  be  used  to  predict  results  for  the  composite,  late  time,  or 
multi-wave  error  and  uncertainty  analysis  based  on  a  multi-path  scattering  formula, 
as  developed  in  [4].  The  input  to  the  multi-path  scattering  error  prediction  formula 
is  elementary,  in  that  it  consists  of  errors  associated  with  individual  isolated  wave 
interactions  and  input  error  uncertainty  only. 

2.1.  The  Single  Propagating  Wave.  We  start  with  the  analysis  of  the 
single  propagating  inward  shock  in  spherical  geometry.  The  radially  dependent 
strength  of  convergent  spherical  shocks  is  discussed  in  [8].  The  spatial  variation  of 
convergent  spherical  shocks  is  contained  in  the  Guderley  solution  [7].  The  inward 
moving  shocks  are  not  of  constant  strength  as  in  a  planar  geometry,  but  evolve  ap¬ 
proximately  according  to  a  power  law  as  a  function  of  the  radius.  From  Whitham’s 
approximation  approach,  we  have 

(2.1)  Mocr"1/n, 
for  cylindrical  shocks,  and 

(2.2)  M  oc  r"2/n, 

for  spherical  shocks.  Here  M  is  the  Mach  number  of  the  shock,  n  =  1  +  ~  + 

and  7  is  the  adiabatic  exponent  defined  as  the  ratio  of  two  specific  heats.  A 
comparison  with  the  exponents  from  Guderley’s  exact  similarity  solution  is  given 
in  Table  2.1.  We  also  have  a  similar  approximate  power  law  for  the  shock  velocity. 

Table  1.  Comparison  of  the  exponents  from  the  approximate  and 
the  exact  similarity  solutions  for  an  inward  propagating  spherical 
shock  wave. 


Cylindrical 

Spherical 

7 

Approximate 

Exact 

Approximate 

Exact 

6/5 

0.163112 

0.161220 

0.326223 

0.320752 

7/5 

0.197070 

0.197294 

0.394142 

0.394364 

5/3 

0.225425 

0.226054 

0.450850 

0.452692 

Fig.  1  shows  the  exponential  divergence  of  the  shock  strength  (here  charac¬ 
terized  by  the  Mach  number)  at  r  -*  0.  The  accuracy  is  amazing  in  view  of  the 
simplicity  of  the  approximate  theory.  The  figure  shows  that  converging  shocks  are 
reacting  primarily  with  the  geometry,  as  assumed  in  the  approximate  theory,  and 
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Figure  1.  Left.  Mach  number  vs.  radius  for  a  single  inward 
propagating  shock.  Right.  The  same  data  plotted  on  a  log- log 
scale. 


are  affected  very  little  by  further  disturbances  from  the  source  of  the  motion;  the 
strength  of  the  initial  shock  enters  only  through  the  constants  of  proportionality  in 
(2.1)  and  (2.2).  This  is  not  true  for  outward  moving  shocks.  They  slow  down  due 
to  both  the  expanding  geometry  and  to  the  continuing  interaction  with  the  flow  be¬ 
hind.  From  Fig.  2,  however,  we  find  that  the  strength  of  an  outward  moving  shock 
also  follows  a  power  law  which  is  similar  to  (2.2)  but  with  a  modified  exponent, 
after  the  radius  of  the  outward  moving  shock  is  three  times  the  initial  radius.  To 
develop  a  model  for  shock  wave  propagation  which  has  a  smaller  pre- asymptotic 
regime,  we  allow  two  distinct  exponents, 


(2.3) 


ro  <r  <  3ro 
r  >  3ro. 


Here  we  choose  ax  =  -0.4,  a2  =  -1.0  for  7  =  1.67,  and  r0  is  the  initial  shock 
radius. 

We  also  study  the  single  propagating  contact  (step  up  and  step  down  cases). 
Fig.  3  shows  the  contact  width  wc  ~  cct x/5  growing  from  2  to  5  cells  with  a  rate 
asymptotically  proportional  to  t1/5.  We  found  that  the  step  up  contact  and  the 
step  down  contact  have  the  same  behavior.  These  properties  appear  to  be  sensitve 
to  the  details  of  the  numerical  algorithm.  We  have  used  a  MUSCL  algorithm.  The 
degree  of  recalibration  of  the  model  as  presented  here,  needed  for  other  solvers, 
is  an  important  question,  out  of  the  scope  of  the  study.  For  some  aspects  of  the 
solution  error,  the  probabilistic  error  formalism  is  more  general  than  is  required. 
When  the  standard  deviation  of  the  error  is  much  smaller  than  the  mean  error 
(when  the  coefficient  of  variation,  their  ratio,  is  close  to  zero),  then  the  error  is 
essentially  deterministic,  and  the  probabilistic  formulation  is  unnecessary.  This  is 
the  case  in  Fig.  3,  with  the  standard  deviation  of  the  width,  shown  to  the  left  scale 
of  Fig.  3,  significantly  smaller  than  the  mean  width. 


2,2.  The  Shock  Contact  Interaction.  We  study  the  wave  strength,  speed, 
width  and  position  errors  after  a  wave  interaction.  We  represent  the  wave  properties 
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Figure  2.  Left.  Mach  number  vs.  radius  for  an  outward  moving 
shock  wave  starting  at  different  radii  rQ.  Right.  The  same  data 
plotted  on  a  log-log  scale;  the  dashed  lines  in  this  plot  represent 
the  power  law  model  (2.3). 


Figure  3.  Ensemble  mean  contact  width  for  a  single  propagat¬ 
ing  contact.  We  record  the  width  in  units  of  Ax.  The  standard 
deviation  is  also  plotted,  as  the  points  to  the  extreme  left. 


as  a  five  tuple 

(2.4)  wak  =  (wak,5l\%,sak,pak), 

where  is  a  wave  strength,  6  is  an  error  in  the  wave  strength,  A  is  a  wave  width,  s 
is  a  wave  speed  error,  and  p  is  a  position  error.  Also  a  =  i  denotes  input  and  a  =  o 
signifies  output.  We  choose  dimensionless  variables  to  measure  wave  strengths;  the 
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Atwood  number  A  =  (p2  —  pi)/(p2  +  p\)  to  measure  the  contact  strength,  and  the 
Mach  number  M  defined  as  the  ratio  of  the  shock  speed  to  the  ahead  state  sound 
speed,  in  the  frame  of  a  stationary  ahead  state,  for  the  shocks.  We  measure  the 
wave  widths  in  units  of  mesh  spacing.  The  wave  position  errors  are  specified  in  grid 
units,  after  an  initial  transient  period.  Within  this  formulation,  we  can  describe 
the  output  wave  errors  by  an  expression  multi-linear  (and  eventually  linear)  in  the 
two  input  wave  strengths,  i.e.  linear  in  each  of  the  two  input  wave  strengths. 

For  input  wave  strengths  u;j,  u/j  (i  denotes  in)  and  output  wave  strengths 
(jj>2  and  o; 3  (o  denotes  out)  (ordered  from  left  to  right),  the  multinomial  expansion 
for  the  output  is  defined  by  its  coefficients  afc.j,  for  J  a  multi  index,  J  —  (j i,  j2)- 
The  expansion  has  the  form 

(2.5)  , 

J 

where  us*'3  =  (u;j)J1(u4)j2’  The  coefficients  depend  parametrically  on  the 
base  case  Riemann  problem,  about  which  a  specified  variation  is  allowed.  Given  a 
statistical  ensemble  of  input  and  output  values  lj1  and  w°,  we  use  a  least  squares 
algorithm  to  determine  the  best  fitting  model  parameters  ak,j ,  for  any  given  poly¬ 
nomial  order  of  model.  We  use  (2.5)  variationally,  that  is  to  map  input  variation 
(about  the  base  case  for  the  ensemble)  to  output  variability.  In  other  words,  (2.5), 
which  is  a  formula  for  wave  strengths,  implies  a  similar  formula  with  different  but 
computable  coefficients  afc,j,  in  which  all  u/’s  are  defined  as  variations  from  the 
base  case,  so  that  they  represent  uncertainty  or  error.  In  the  planar  case  [4],  we 
showed  that  a  linear  input-output  relation,  —  ctk,o  +  was  sufficient  to 

describe  the  exact  (nonlinear)  input-output  relation.  The  adequacy  of  this  model 
is  established  from  consideration  of  the  STD  of  the  model  error.  We  see  that  the 
same  assumption  is  adequate  in  the  present  case.  We  also  have  an  expansion  similar 
to  (2.5)  for  the  wave  strength  errors,  wave  width  errors  and  wave  position  errors. 
To  avoid  possible  confusion,  we  note  that  the  word  error  is  used  to  describe  two 
different  types  of  quantities.  The  solution  (wave  strength,  width  or  position)  error 
is  a  difference  between  a  fine  grid  and  a  coarse  grid  solution.  The  coefficients  for 
a  linear  model  of  point  values  of  this  error  are  presented.  The  model  error  is  the 
difference  in  some  quantity  and  the  linear  model  (2.5)  to  represent  it,  as  a  function 
of  the  variable  parameters  which  define  the  ensemble. 

We  begin  with  the  analysis  of  the  initial  shock  contact  SNRP  at  the  ensemble 
averaged  level.  We  present  the  linear  model  coefficients  in  Table  2,  with  ±10% 
variation  for  the  initial  contact  strength  and  ±5%  variation  for  the  initial  shock 
strength  (consistent  with  ±10%  variation  in  pressure  ratio  as  used  in  the  planar 
study  [4])  about  the  base  case.  According  to  the  analysis  of  Sec.  2.1,  the  strength 
of  this  initial  inward  shock  is  not  constant,  and  is  increasing  as  it  moves  toward  the 
origin.  We  use  the  power  law  M  —  Cr~2/n  to  estimate  the  initial  shock  strength 
at  the  interaction  time  and  use  this  quantity  represented  by  the  variable  C  as  the 
input  shock  strength  in  the  modeling.  The  input  contact  width  has  been  set  to 
zero,  as  part  of  the  specification  of  this  SNRP. 

To  read  Table  2,  we  note  that  the  first  (a;")  row  (labeled  in  the  table  as 
wf  (1.  sonic))  lists  coefficients  for  J  —  (0,0),  J  —  (1,0),  etc .  These  coeffi¬ 
cients  are  determined  by  a  least  squares  algorithm  that  minimizes  the  expected,  or 
mean  error  over  the  ensemble,  in  comparing  the  linear  predictions  to  the  exact  so¬ 
lution  of  the  Riemann  problem.  The  last  two  columns  describe  errors  in  the  model 
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Table  2.  The  SNRP  shock  contact  interaction.  Expansion  coeffi¬ 
cients  for  output  wave  strengths,  wave  strength  errors,  wave  width 
errors  and  wave  position  errors  (linear  model)  for  the  initial  shock 
contact  interaction.  Here  the  base  case  input  contact  wave  width 
is  zero.  The  final  columns  refer  to  difference  between  the  linear 
model  (2.5)  and  the  exact  quantity.  The  errors  in  rows  4-12  refer 
to  the  difference  between  the  numerical  solution  on  100  cells  and 
the  exact  solution  using  2000  cells. 


variable  \  coef 

const 

Ul{ 

U>2 

model  error 

(contact) 

(1.  sonic) 

STD 

STD/J° 

wave 

strengths  (100  cells) 

u>°  (1.  sonic) 

-33.353 

19.521 

2.501 

0.860 

0.954% 

u>2  (contact) 

0.374 

0.200 

0.0003 

0.042 

7.650% 

u>3  (r.  sonic) 

3.568 

0.402 

-0.045 

0.009 

0.463% 

wave  strength  errors  (100  cells) 

S i  (1.  sonic) 

2.039 

-3.200 

-0.01 

0.157 

0.174% 

#2  (contact) 

0.236 

0.016 

-0.002 

0.021 

3.825% 

S\ ij  (r.  sonic) 

0.053 

0.003 

-0.001 

0.0008 

0.041% 

wave  width  errors  (100  cells) 

Xi  (1.  sonic) 

1.675 

0.305 

0.017 

0.085 

A£  (contact) 

7.093 

0.482 

-0.146 

0.239 

(r.  sonic) 

2.829 

0.302 

-0.024 

0.107 

wave  position  errors  (100  cells) 

p °  (1.  sonic) 

-0.247 

0.242 

0.005 

0.009 

P2  (contact) 

0.643 

0.065 

-0.011 

0.192 

(r.  sonic) 

-0.042 

0.062 

0.004 

0.009 

(2.5).  The  presence  of  outliers  was  monitored  and  the  ensemble  L norm  deter¬ 
mined  (results  not  tabulated);  occasional  outliers  indicate  non  Gaussian  statistics. 
The  model  error  is  defined  as  (predicted  -  exact)  where  exact  is  the  result  of  the 
simulation  and  predicted  is  the  value  given  by  the  finite  polynomial  (linear)  model 
(2.5).  The  column  STD  is  the  standard  deviation  of  (predicted  -  exact).  Note 
that  the  STD  errors,  as  defined,  are  dimensionful.  To  aid  in  interpreting  the  error 
magnitudes,  we  present  in  a  final  column  (labeled  STD/cu°)  the  standard  deviation 
of  the  error  in  the  model  divided  by  the  mean  value  of  the  variable  predicted.  This 
column  represents  a  fractional  (dimensionless)  error  in  the  model. 

According  to  the  analysis  of  Sec.  2.1,  the  strength  of  the  output  inward  moving 
shock  is  modeled  as  CV~2/n.  This  formula  is  accurate  after  some  time,  and  the 
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Table  2  entry  is  w°  =  C  in  this  formula.  We  form  a  linear  model  for  this  constant 
in  this  expression  in  Table  2.  We  find  very  small  errors  in  the  exponent,  not 
tabulated  here.  We  developed  a  model  (2.3)  for  the  strength  of  the  output  outward 
moving  shock  in  Sec.  2.1.  Here  in  our  study,  we  are  only  concerned  with  the  first 
formula  in  (2.3).  The  entry  u)%  in  the  table  represents  the  coefficient  multiplying 
the  power  term. 


Figure  4.  Left:  ensemble  mean  inward /outward  -moving  shock 
and  contact  widths  after  a  shock  contact  interaction.  Right:  en¬ 
semble  mean  shock  and  contact  position  errors  as  a  function  of 
time,  expressed  in  grid  units.  The  associated  standard  deviations 
are  extremely  small,  not  shown  in  the  plots.  In  the  legend,  C.  de¬ 
notes  the  contact  while  I.S.  and  O.S.  are  the  inward  and  outward 
moving  shocks. 

The  three  variable  (A)  rows  in  Table  2  represent  wave  width  errors.  The  stan¬ 
dard  deviation  for  this  quantity  is  about  10%  of  the  mean  value,  indicating  that 
the  error  model  is  (on  the  whole)  satisfactory,  and  that  the  shock  wave  widths  are 
not  (mostly)  fluctuating  greatly.  The  inward  moving  shock  width  decreased  about 
10%  relative  to  the  wave  width  at  the  interaction  time,  while  the  outward  moving 
shock  width  increased  about  10%.  See  Fig.  4,  left  frame.  The  contact  width  is 
modeled  as  cct x/5  where  both  the  width  and  t  are  expressed  in  mesh  units.  The 
Table  2  entry  A£  =  cc  in  this  formula.  We  form  a  linear  model  for  this  constant  in 
this  expression  in  Table  2. 

We  also  study  the  wave  position  errors.  Fig.  4,  right  frame,  shows  the  position 
errors  as  a  function  of  time.  The  entries  in  the  wave  position  rows  of  Table  2 
present  those  errors,  given  in  mesh  units.  All  position  errors  are  subgrid.  The 
standard  deviations  are  smaller  than  the  means,  indicating  that  the  errors  are 
basically  deterministic. 

All  solution  errors  are  sensitive  to  the  grid  spacing,  taken  to  be  100  computa¬ 
tional  cells  in  Table  2-4.  This  sensitivity  is  not  extreme.  For  example,  if  the  100  cell 
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model  is  used  to  analyze  the  500  cell  data,  the  model  errors  (STD)  approximately 
double,  but  remain  small. 

2.3.  Shock  Reflection  at  the  Origin.  Here  we  study  the  reflection  of  the 
shock  off  the  origin.  According  to  the  analysis  of  Sec.  2.1,  the  input  inward  moving 
shock  has  infinite  strength  at  the  origin.  We  used  the  strength  at  the  radius  r  =  1 
as  the  initial  state  and  the  input  shock  strength  in  the  modeling  process.  We  study 
the  wave  strength,  wave  strength  errors,  wave  width  and  wave  position  errors.  See 
Table  3. 


Table  3.  The  SNRP  defined  by  the  shock  reflection  at  the  origin. 
Expansion  coefficients  for  output  wave  strengths,  wave  strength 
errors,  wave  width  errors  and  wave  position  errors  (linear  model) 
for  input  variation  ±10%. 


variable  \  coef 

const 

u>\ 

(1.  sonic) 

model  error 

STD  STD/J0 

u°  (r.  sonic) 

-242.394 

5.606 

1.137 

0.468% 

5°  (r.  sonic) 

-3.27 

0.031 

0.112 

0.045% 

A°  (r.  sonic) 

1.221 

0.018 

0.099 

Pi  (r.  sonic) 

0.474 

0.001 

0.012 

We  found  that  the  Mach  number  of  the  outward  moving  shock  (reflected  shock) 
was  essentially  independent  of  the  input  variation  in  Mach  number.  To  explain  this 
phenomena,  we  recall  that  the  ambient  state  ahead  of  the  outward  moving  reflected 
shock  is  an  incoming  continuously  variable  flow.  The  sound  speed  ahead  of  this  flow 
is  affine  linearly  dependent  on  the  strength  of  the  incoming  shock  wave,  as  is  the 
shock  speed  of  the  reflected  outward  moving  shock  wave.  Thus  the  outward  moving 
Mach  number,  as  a  ratio  of  two  quantities  varying  affine  linearly  with  the  incoming 
shock  strength,  has  a  fractional  linear  form  in  the  incoming  wave  strength.  A  simple 
calculation  shows  that  the  variation  in  the  outward  moving  shock  Mach  number 
M0  contains  the  factor  (1  -  Ma)  and  since  M0  «  1.2,  this  small  factor  suppressed 
variation  in  M0  as  a  function  of  Af*,  the  Mach  number  of  the  incoming  shock.  Thus 
the  Mach  number  is  not  a  good  measure  for  the  outward  moving  shock  strength. 
We  choose  the  pressure  behind  the  reflected  shock  instead  as  u°  in  Table  3.  The 
pressure  also  follows  the  power  law.  The  large  entries  in  this  row  result  from  the 
fact  that  the  (dimensional)  pressure  (u>?)  is  much  larger  in  pressure  units  than  the 
Mach  number 

2.4.  The  Contact  Reshock  Interaction.  After  reflection  from  the  origin, 
the  transmitted  lead  shock  wave  re-crosses  the  deflected  contact.  The  outgoing 
waves  from  this  interaction  consist  of  a  rarefaction  wave  propagating  toward  the 
origin,  a  contact  and  a  shock  propagating  outward.  The  region  inside  of  the  outward 
propagating  shock,  on  both  sides  of  the  contact  is  not  piecewise  constant,  but 
contains  an  inward  propagating  compression,  which  eventually  breaks  to  form  an 
inward  moving  shock,  reaching  the  origin  at  interaction  4.  This  inward  moving 
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Table  4.  The  SNRP  contact  reshock  interaction.  Expansion  co¬ 
efficients  for  output  wave  strengths,  wave  strength  errors,  wave 
width  errors  and  wave  position  errors  (linear  model). 


variable  \  coef 

const 

u4 

U>2 

model  error 

(r.  sonic) 

(contact) 

STD 

STD/u>° 

wave 

strengths  (100  cells) 

u°  (1.  sonic) 

0.097 

-0.108 

0.436 

0.031 

13.305% 

u>2  (contact) 

0.103 

-0.192 

1.168 

0.007 

1.116% 

u>3  (r.  sonic) 

0.988 

0.195 

-0.225 

0.003 

0.262% 

wave  strength  errors  (100  cells) 

(1.  sonic) 

-0.291 

0.161 

-0.468 

0.017 

7.296%  1 

82  (contact) 

-0.067 

0.142 

-0.125 

0.006 

0.957% 

<$3  (r.  sonic) 

-0.030 

0.107 

-0.0004 

0.001 

0.087% 

wave  width  errors  (100  cells) 

A?  (1.  sonic) 

9.776 

-6.372 

5.091 

0.484 

Ag  (contact) 

1.903 

0.156 

-0.677 

0.534 

Ag  (r.  sonic) 

4.088 

-1.401 

1.549 

0.168 

wave  position  errors  (100  cells) 

p°  (1.  sonic) 

4.782 

-3.602 

2.372 

0.379 

P2  (contact) 

-0.453 

0.409 

-0.054 

0.177 

pg  (r.  sonic) 

-0.199 

-0.685 

3.213 

0.052 

compression  is  generated  from  the  geometrically  caused  weakening  of  the  outward 
moving  shock,  and  is  a  well  recognized  aspect  of  spherical  shock  wave  dynamics. 
The  shock  and  the  rarefaction  interact,  and  eventually  the  rarefaction  disappears  in 
this  interaction.  Here  we  only  follow  the  waves  through  the  output  of  interaction  3, 
and  thus  avoid  much  of  this  interaction.  Specifically,  we  focus  on  the  inward  moving 
rarefaction  and  not  the  inward  moving  shock.  We  study  the  wave  strength,  wave 
strength  errors,  wave  width  and  wave  position  errors  resulting  from  interaction  3. 
See  Table  4. 

According  to  the  analysis  of  Sec.  2.1,  this  is  a  step  down  interaction  and  the 
contact  width  is  modeled  as  cct1^  where  both  the  width  and  t  are  expressed  in 
mesh  units.  We  form  a  linear  model  for  the  coefficient  cc  in  this  expression  in 
Table  4.  The  rarefaction  width  has  the  form  constant  +  rate  x  time.  The  entry 
A?  refers  to  the  constant,  which  gives  an  offset  for  the  centering  of  the  rarefaction 
wave.  This  entry  is  expressed  in  mesh  units. 
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3.  Composite  Shock  Interaction  Problems 


i _ i _ i _ i _ I — i — L 


Figure  5.  Left:  space  time  density  contour  plot  for  the  multiple 
wave  interaction  problem  studied  in  this  section,  in  spherical  ge¬ 
ometry.  Right:  type  and  location  of  waves  determined  by  the  wave 
filter  analysis  with  labels  for  the  interactions.  Here  I.R.  denotes 
the  inward  moving  rarefaction. 

We  consider  the  repeated  interactions  of  a  spherically  symmetric  shock  wave 
with  a  spherical  contact  located  near  the  origin.  The  base  case  for  each  wave 
interaction  coincides  with  the  base  case  assumed  for  the  interactions  studied  in 
Sec.  2.  The  transmitted  shock,  after  interaction  with  the  contact,  progresses  to 
interact  with  ( i.e .  reflect  off)  the  origin.  This  interaction  was  also  studied  in 
Sec.  2.  Subsequently,  there  are  a  number  of  reverberations,  of  reflected  rarefactions 
and  compression  waves,  between  the  contact  and  the  origin.  The  interactions  are 
illustrated  by  the  space  time  contour  plots  of  the  density,  shown  in  Fig.  5,  left.  In 
Fig.  5  (right),  we  show  the  type  and  location  of  the  waves,  as  determined  by  the 
wave  filter  analysis  program.  Five  interactions  are  extracted  from  the  complex  wave 
interaction  problem.  Both  figures  refer  to  the  base  case.  The  build  up  of  complex 
wave  patterns  is  evident. 

The  main  point  of  this  section  is  to  formulate  and  validate  the  multipath  scat¬ 
tering  formula  of  [4]  for  analysis  of  errors.  We  analyze  errors  at  the  output  to 
interaction  3  directly,  comparing  the  100  mesh  and  500  mesh  simulation  to  a  2000 
mesh,  fine  grid  simulation,  here  taken  as  a  substitute  for  the  exact  solution.  These 
errors  are  compared  to  those  generated  by  adding  up  and  propagating  errors  from 
the  input  data  and  from  each  of  the  interactions  1  to  3,  using  the  multipath  scat¬ 
tering  formula.  Thus,  for  example,  a  position  error  as  input  to  interaction  1  is 
translated  geometrically  to  a  position  error  for  the  output  to  interaction  1  via  sim¬ 
ple  geometric  considerations  as  in  [4] .  This  error  is  propagated  to  an  input  error  for 
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Table  5.  Predicted  and  simulated  errors  for  output  wave 
strengths,  wave  widths  and  wave  positions,  output  to  interaction 
3.  The  inward  rarefaction  and  contact  strengths  are  expressed  di- 
mensionlessly  as  Atwood  numbers.  The  outward  shock  strengths 
are  in  the  units  of  Mach  number.  The  width  and  position  errors 
are  in  mesh  units.  The  wave  strength  errors  are  expressed  as  mean 
±  2  a  where  a  is  the  ensemble  STD  of  the  err  or /uncertainty. 


variable  \  error 

Simulation  Prediction 

100  vs.  2000  mesh 

Simulation  Prediction 

500  vs.  2000  mesh 

wave  strength  errors  and  propagated  initial  uncertainties 

5?  (1.  sonic) 

0.04±2(0.03) 

0.03±2(0.02) 

0.01±2(0.02) 

0.009±2(0.01) 

<$2  (contact) 

0.14±2(0.05) 

0.12±2(0.02) 

0.03±2(0.01) 

0.03±2(0.008) 

Si jj  (r.  sonic) 

-0.02±2(0.02) 

-0.02±2(0.01) 

-0.006±2(0.005) 

-0.007±2(0.004) 

mean  wave 

width  errors 

mean  wave 

width  errors 

A?  (1.  sonic) 

3.04 

2.83 

2.63 

2.72 

A£  (contact) 

5.36 

6.11 

5.56 

6.08 

A§  (r.  sonic) 

2.71 

3.04 

2.92 

2.98 

mean  wave  position  errors 

mean  wave  position  errors 

(1.  sonic) 

1.25 

0.23 

0.12 

0.18 

P2  (contact) 

0.43 

0.06 

0.05 

0.04 

p3  (r.  sonic) 

-0.73 

-0.15 

-0.08 

-0.11 

interaction  2  through  solutions  of  radial  differential  equations.  Propagation  con¬ 
tinues,  and  yields  an  error  at  the  output  to  interaction  3.  See  Table  5.  The  wave 
strength  rows  present  the  result  of  initial  uncertainty  propagated  to  the  output  of 
interaction  3  as  well  as  the  accumulation  of  solution  errors.  The  multipath  scat¬ 
tering  formula  gives  reasonable  prediction  of  error  magnitudes  in  all  cases  except 
the  wave  position  errors  for  the  under  resolved  (100  mesh)  simulation.  We  see  that 
the  created  numerical  solution  errors  are  important.  We  also  find  that  a  major 
portion  of  the  created  numerical  solution  errors  come  from  the  second  interaction, 
the  shock  reflection  interaction.  A  detailed  study  of  these  errors  and  their  relative 
importance  will  be  presented  in  a  following  paper. 

4.  High  Dimensional  Wave  Filter 

4,1.  The  Two  Pass  Algorithm.  We  develop  a  two  dimensional  wave  filter. 
The  analysis  starts  at  a  general  point,  and  at  that  point,  searches  in  a  general 
direction.  The  analysis  proceeds  in  two  passes.  The  first  pass  uses  only  state  data 
along  a  line  segment,  which  are  extracted  from  the  2D  simulation  data  and  analyzed 
using  the  ID  filter  introduced  previously  [4].  In  this  first  pass,  we  choose  the  normal 
direction  as  the  direction  through  the  given  point  which  has  the  strongest  wave  in 
the  dominant  family  and  the  weakest  waves  in  the  other  families.  We  obtain  an 
approximate  wave  position  and  width  by  applying  the  ID  filter  to  the  line  segment 
in  the  normal  direction  of  the  wave.  The  second  pass  improves  in  this  choice  of 
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Figure  6.  The  strengths  of  3  output  waves. 

normal  direction  by  the  construction  of  a  local  tangent  plane.  We  introduce  an  array 
of  (two  in  2D)  additional  ID  line  segments  displaced  tangentially  ( i.e .  normal  to 
the  normal  just  constructed)  from  the  original  one.  The  analysis  as  above  along 
each  segment  in  the  array  leads  to  wave  positions  in  the  tangentially  displaced 
array.  Finally  a  tangent  plane  is  fit  to  array  of  wave  positions.  The  normal  to  this 
plane  is  the  wave  normal  of  the  final  construction.  We  analyze  the  inward  moving 
shock  wave  before  interaction  1  of  our  base  case  using  this  wave  filter. 

4.2.  Initial  Estimate  for  the  Normal  Direction  of  the  Wave.  Applying 
a  ID  wave  filter  to  the  state  data  along  a  line  segment  gives  a  wave  position  and  the 
width.  Starting  from  a  trial  wave  position,  and  a  trial  normal  direction  given  by 
an  angle  9 ,  and  a  distance  r  along  the  resulting  segment,  we  pick  up  two  points  on 
the  left  and  the  right  ends  of  the  segment.  With  the  suitable  choice  of  r,  (e.g.  the 
wave  width)  the  states  at  the  two  points  represent  the  left  and  right  constant  states 
of  the  wave.  The  left  and  right  states  and  their  decomposition  into  ID  elementary 
wave  via  a  Riemann  solution  varies  according  to  9.  The  Riemann  solution  has 
three  output  waves,  whose  strengths  are  assessed  dimensionlessly.  The  type  of 
the  strongest  wave  is  defined  to  be  the  wave  type  at  the  trial  wave  position.  The 
strongest  wave  has  a  maximum  value  and  the  weaker  waves  have  minimum  value 
at  an  angle  $m.  We  choose  9m  as  the  normal  direction.  Using  this  normal  direction 
9m >  the  ID  wave  filter  will  select  a  trial  position  pm  and  a  wave  width.  Fig.  6  shows 
the  strengths  of  three  output  waves  at  an  angle  9, 0  <  9  <  90°. 

4.3.  Corrected  Estimate  for  the  Normal  Direction  of  the  Wave.  The 
accuracy  of  the  wave  position  and  normal  is  improved  through  the  construction  of 
a  local  patch  of  a  wave  surface.  The  2D  wave  surface  patch  is  a  curve  segment 
defined  by  two  linear  segments  meeting  at  Pm,  and  thus  defined  (in  the  present 
approximation)  by  two  points  Pi  and  Pr ,  serving  as  end  points  of  the  two  linear 
segments.  We  construct  trial  points  P/,  P/  by  their  locations  along  the  tangential 
directions  9m  -  90 °,9rn  +  90°  respectively  passing  through  Pm  at  a  certain  distance 
d.  Two  normal  direction  9i,9r  at  P/,P/  respectively  are  found  as  above.  Using 
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Figure  7.  Schematic  diagram  illustrating  the  computation  of  the 
corrected  estimate  for  the  normal  to  the  wave. 

the  ID  (first  pass)  construction  along  the  lines  defined  by  Oi ,  P/  and  6r,  P/,  we 
find  two  points  Pi,Pr  on  the  interface  to  define  the  local  wave  surface  patch.  With 
these  three  points,  Pi,  Pm,  Pr  on  the  interface,  we  can  construct  a  circle.  This  circle 
represents  the  interface  locally.  The  normal  direction  0n  of  the  circle  at  the  point 
Pm  represents  the  normal  direction  of  the  interface  at  the  point  Pm.  We  choose  0n 
as  the  normal  direction  for  the  correction  step. 

Fig.  7  shows  a  schematic  diagram  illustrating  the  correction  step.  We  see  the 
correction  from  the  initial  estimate  lni  to  the  corrected  normal. 

The  correction  depends  on  the  distance  d  from  Pm  to  P/  and  P/.  If  d  is  too 
small,  the  normal  direction  is  sensitive  to  the  error  in  the  points  Pi,Pr  on  the 
interface.  Fig.  8  shows  that  the  choice  of  d  equal  to  3  mesh  units  is  satisfactory. 
This  figure  also  illustrates  the  improvement  of  accuracy  and  the  robustness  of  the 
correction  step.  Here,  we  tested  89  lines  at  angles  1  —  89°  for  Oq  from  (0,0)  for  Po 
to  the  initial  shock  data,  input  to  interaction  1  of  our  base  case  simulation.  We 
see  that  the  standard  deviations  of  the  errors  for  the  normal  direction  are  small. 
The  space  averages  of  the  errors  for  the  normal  direction,  which  we  do  not  present, 
were  very  small  (less  than  10%  of  the  STD)  for  all  choices  of  d .  All  choices  of  d 
give  a  satisfactory  normal  direction.  The  initial  errors  are  somewhat  large.  This 
is  because  the  wave  has  0  wave  width  at  the  initial  time  step.  The  ID  filter  fits 
state  data  to  an  error  function,  and  this  fit  is  not  stable  before  the  initialization 
transients  in  the  numerical  shock  have  disappeared. 

5.  Conclusion 

The  multipath  integral  formula  to  express  solution  error  as  a  composition  of 
errors  associated  with  individual  interactions  and  a  propagation  law  for  errors  be¬ 
tween  interactions  has  been  extended  to  spherically  symmetric  shock  interaction 
problems.  The  main  new  difficulties  encountered  were  the  non-constancy  of  the 
solution  between  interaction  events  and  the  non-constancy  of  waves  and  errors  be¬ 
tween  interactions.  In  addition,  as  preparation  for  future  multidimensional  error 
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Figure  8.  STD  of  errors  for  the  wave  normal  constructed  from 
the  2D  wave  filter,  as  a  function  of  the  tangential  separation  d  of 
the  wave  front  points  P[yP^  in  mesh  units. 


analysis,  we  have  extended  the  wave  filter  previously  introduced  to  detect  waves  in 
multidimensional  numerical  solutions. 
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